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A method for arbitrary two-dimensional scaling of phonocardiography signals. 



5 The invention relates to a method for assisting medically trained persons in 

diagnosing and classifying medical conditions of the heart by means of auscultation. 
It also relates to a method and apparatus for transforming the information content of 
a phonocardiography signal to a frequency range and duration commensurate with a 
communication channel having properties that fall outside the normal frequency 
10 range of phonocardiography signals. 



Signals obtained by means of a transducer are phonocardiography representations of 
sounds traditionally listened to by means of a stethoscope. Training in auscultation 
takes a long time and requires an aptitude for recognising and classifying aural cues, 
1 5 frequently in a noisy environment. 20-30 different conditions may need to be 

differentiated, and within each, the severity evaluated. Furthermore, there may be 
combinations among these. These factors contribute to explaining why not all 
physicians perform equally well when diagnosing heart conditions, and why it may 
be time-consuming. 

20 

Use of modern signal processing has been proposed to enhance the relevant aural 
cues, and in some cases it has been demonstrated that a transposition into a different 
frequency range or repetition of a repetitive signal at a different repetition frequency 
from that originally recorded may be advantageous. Among the patents that describe 
25 this kind of signal processing the following may be mentioned: 



US 4,220,160 uses the heart sound picked up by a microphone to modulate a higher 
and more audible carrier frequency. US 4,424,815 is a very fundamental patent 
employing a range of solutions, from magnetic drum (i.e. a non-digital procedure) to 
30 a tapped delay line for stretching time while conserving frequency content. US 

4,528,689 describes a time stretching digital stethoscope that records several cycles, 
and obtains a "slow-down factor" proportional to the number of cycles recorded. 
Each part of a sound between zero crossings is repeated with opposite polarity, 
thereby obtaining a signal that sounds like original signals but lasting longer. Only 
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whole-number factors are possible. US 4,792,145 describes how frequency scaling is 
5 obtained by multiplying each Fourier spectral component by a factor. The original 
relationship between frequency components is retained, and it is purported that the 
audibility is improved. Similarly, the time may be scaled. 

All of the above solutions have severe limitations, either because of an increase in 
1 0 non-linear distortion and noise, or because the factor for stretching can only be 
whole-numbered. 



It has been proposed to use an algorithm based on Matching Pursuit, such as in Xuan 
Zhang et al. "Time-Frequency scaling transformation of the phonocardiogram based 

15 on the matching pursuit method", IEEE Trans. Biom. Eng. Vol. 45, No. 8, pp. 972- 
979 in view o/Xuan Zhang et al. "Analysis-synthesis of the phonocardiogram based 
on the matching pursuit method", idem pp. 962-971 (August 1998). However, this 
procedure is an iterative algorithm that relies on the modulation of the signal as a 
sum of functions selected from a library of functions that has to be searched several 

20 times for each signal segment. It is hence a very time-consuming procedure, and in 
practice the control over the individual phases of the segments required for the 
stability of the scaled signal is questionnable. 

According to the invention, these disadvantages may be avoided by synthesizing a 
25 "clean" signal consisting of a sum of sinusoids and to perform a transformation on 
these sinusoids in order that the relative amplitudes are maintained while the 
frequency range or time axis are changed in a continuous range of scales. In a further 
refinement of this method, the scale may be linked to other phenomena by means of 
an auto-scale function. 

30 

An embodiment is particular in that that the signals are converted on a running basis 
to a sinusoidal model representation, that a time and/or frequency axis scaling is 
defined and used to control the amplitudes and phases of the sinusoids, which are 
subsequently added to create a time and/or frequency scaled representation of the 



phonocardiography signal. The basis of the invention is the method in which a 
correct amplitude and phase adjustment of the contributing sine generators is 
obtained. 

5 An advantageous embodiment of the invention is particular in that it comprises the 
steps of obtaining the frequency content of the signal by applying a Short Time 
Fourier Transform with overlapping segments, a Discrete Fourier Transform being 
performed on each segment, performing a frequency peak search on each segment by 
consecutive removal of the spectral components having the highest energy content, 

10 identifying each peak by its frequency value repeating the peak search, until a 

maximum number of peaks have been identified or until the energy content of the 
last peak is below a preset minimum, establishing a segment-by-segment map of 
spectral peaks, said peaks forming a track over time, optionally subjecting the 
frequency values of the spectral peaks to a multiplication, defining a synthesis frame 

15 of time, based on said segments, optionally subjecting each frame to a multiplication 
of the time scale, adjusting the phase of sine generators centered on the frequencies 
of the tracks, adjusting the amplitudes of said sine generators, and summing the 
outputs of all sine generators active at any one instant for a given frame length T. 

20 A more general embodiment is particular in that it includes the following steps: 
obtaining the frequency content of the signal by applying a Short Time Fourier 
Transform with overlapping segments, a Discrete Fourier Transform being 
performed on each segment, performing a frequency peak search on each segment by 
repeated removal of the highest spectral "hills" identifying each peak by its 

25 frequency value, repeating the peak search until a maximum number of peaks have 
been identified or until the peak level of the last peak is below a preset minimum, 
establishing a segment-by-segment map of spectral peaks, said peaks forming a track 
over time, adjusting the phase of sine generators centered on the frequencies 
represented by the tracks, summing the outputs of all sine generators active at any 

30 one instant for a given frame length T, creating a continuous output signal by joining 
consecutive frames. In the present embodiment use is made of the fact that the 
estimate of the importance of a peak (a "hill") may be obtained by other means than 
an energy measure, i.e. a different sorting criterion. 
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A further advantageous embodiment for scaling a phonocardiography signal on 
the frequency axis is particular in that the spectral peaks are multiplied by a factor q. 

A further advantageous embodiment for scaling a phonocardiography signal on the 
5 time axis by a desired factor is particular in that the frame length is multiplied by a 
factor p. 

A further advantageous embodiment for autoscaling a phonocardiography signal on 
the time axis is particular in that the scaling factor p is set such that the frame length 
10 multiplying factor is equal to the heart rate divided by 60. This means that a stability 
of the signal will be obtained irrespective of changes in the heart rate. 

An advantageous embodiment of the invention that ensures sufficient resolution and 
precision in synthesis is particular in that the number of sine generators is maximum 
15 50 in any one frame. 

The invention also relates to an apparatus for performing the method, being 
particular in comprising means for windowing the time function, short time Fourier 
spectrum analysis means, means for searching and classifying spectral peaks, means 

20 for comparing phases of signals corresponding to said spectral peaks, means for 
controlling the phases of sine generators providing signals corresponding to said 
spectral peaks, means for controlling the amplitudes of said sine generators, and 
means for summing the signals of said controlled sine generators in order to obtain a 
synthesized and essentially noise free output signal representative of said time 

25 function. 



The invention will be more fully described in the following with reference to the 
drawing, in which 

30 Fig. 1 shows a block diagram representative of the operations performed on an input 
signal in order that a sinusoidal model is obtained for the signal, which is 
summed to obtain an output signal, 
Fig. 2 shows the procedure of block segmentation, 
Fig. 3 shows the estimation of a finite number of spectral peaks, 
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Fig. 4 shows results of the peaks search 
Fig. 5 shows the matching procedure 
Fig. 6 shows the final results of the peak matching 
Fig. 7 shows boundary condition for the interpolation functions, 
5 Fig. 8 shows interpolation of instantaneous phase function and amplitude function 
Fig. 9 shows a spectrogram of the original heart sound, 
Fig. 10 shows a spectrogram of the sinusoidal coded heart sound, 
Fig. 1 1 shows a tracking pattern for the sinusoidal coded heart sound, 
Fig. 12 shows a Time-Frequency scaling procedure, 
10 Fig. 13 shows a spectrogram of a frequency shifted heart sound (p=l,^=l,r=100Hz), 
Fig. 14 shows a spectrogram of a frequency scaled heart sound (p=l,<7=2,r=0), 
Fig. 15 shows a spectrogram of a non-scaled sinusoidal coded heart sound 
(p=l,q=l,r=0), 

Fig. 16 shows a spectrogram of a joint time-frequency scaled heart sound (p=2,<7=2, 
15 r=0), 

Fig. 17 shows a spectrogram of a time scaled heart sound (p=2,q=l,r=0), and 
Fig. 18 shows a spectrogram of the original heart sound (p=l,q=l,r=0), but using a 
different time axis. 

20 The mathematical foundation for this procedure is described in the following 
paragraphs. 

The sinusoidal model was originally developed for speech applications by McAulay 
and Quantieri and a thorough description can be found in Robert J. McAulay and 

25 Thomas F. Quatieri, "Speech Analysis/Synthesis Based On A Sinusoidal Representa- 
tion", IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP- 
34, NO. 4, August 1986, page 744-754, and Thomas F. Quatieri and Robert J. 
McAulay, "Speech Transformations Based On A Sinusoidal Representation", IEEE 
Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-34, NO. 6, 

30 December 1986, page 1449-1464. Because the spectral structure of a speech signal 
and a heart signal are very different (a heart beat does not have pitch like voiced 
speech, nor a broadband structure like unvoiced speech), the original sinusoidal 
model will not be able to model a heart beat properly. A sinusoidal model specific 
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for heart sounds is developed by exploring the spectral characteristics of 

the heart beat, when tracking spectral components in the signal. 

The sinusoidal model is based upon sinus functions. The synthesis part of the model 
5 consists in adding a group of sinusoidal functions: 



Eq.(l): 

y(0 = X^(Ocos(n*(0) 

k - 1 

The number of functions L(t) will change as a function of time, depending on the 
complexity of the signal. Each sinusoidal function (frequency track) is controlled by 
10 the total phase function Q.(t), and the amplitude A(t) . The main objectives of the 
analysis part of the model is to establish the content of these two functions. 

A fundamental relation for the sinusoidal model is the relationship between total 
phase and the instantaneous frequency: 



15 Eq. (2): 



co(0 = T n (0 [rad/s] 
at 

t 

Q(0 = \ (o(x)dx + q> offset [rad] 

—CO 

Substituting the expression for the total phase into eq. (1), the final expression for the 
sinusoidal model is obtained: 



Eq. (3): 

L(t) f t 

y(t) = £^(0cos 



J co k (z)dT + <p k 



V-oo 



20 The structure of a traditional sinusoidal model is illustrated in Fig. 1 . The input 

signal x(t) is obtained from a sensor for heart sounds placed in contact with the body 
under examination. The output signal y(f) is made available to the medically trained 
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person via e.g. headphones after suitable amplification. Furthermore it may be 
used for a visual display and comparison of frequency spectra on a display screen. 

The model may be divided into three parts: an Analysis part where the signal 
parameters are measured, an Interpolation and transformation part, where Q(t) and 
5 A(t) are interpolated from one window position to the next window position, and a 
Synthesis part where the transformed heart signal is reconstructed. 

The apparatus, the principle of which is shown in Fig. 1 , essentially comprises means 
for windowing the time function, short time Fourier spectrum analysis means 1, 

10 means for searching and classifying spectral peaks 2, 3, means 4 for comparing 
phases of signals corresponding to said spectral peaks, means for controlling the 
phases of sine generators 6 providing signals corresponding to said spectral peaks, 
means 5, 7 for controlling the amplitudes of said sine generators, and means 8 for 
summing the signals of said controlled sine generators in order to obtain a 

1 5 synthesized and essentially noise free output signal representative of said time 

function. The content of each block of Fig. 1 will be elaborated more in the following 
sections. 

The first step of the Analysis is to measure the frequency content of the signal. As 
20 the frequency content is changing as a function of time, a Short Time Fourier 

Transform (STFT) is used to track the frequency content. The STFT segments the 
signal into overlapping segments (75%), and each segment is Fourier transformed 
using a Discrete Fourier Transform (DFT). The segmentation procedure is illustrated 
in Fig. 2. 

25 

When a signal segment has been Fourier transformed, the most dominant frequency 
component in the spectrum has to be found, and used as local "check points" for a 
frequency track over time. The original sinusoidal model finds every peak in the 
spectrum, and associates a sinusoidal generator to each peak. In the case of a voiced 
30 speech signal, the signal is periodic and the corresponding spectrum will be discrete 
(pitch). As the spectrum is discrete, the known peak search procedure will find a 
large number of peaks. However, a heart beat is non-periodic and has a nature like a 
Gaussian function modulated by a low- frequency sinusoid. The corresponding 



8 

spectrum is also Gaussian, and the known peak search procedure will model this 
heart beat by only one sinusoidal function, where the amplitude is equal to the peak 
level. This will not model the heart beat correctly, compared to the remaining 
spectrum. However, if the amplitude of the sinusoidal function instead is set equal to 
5 the area below the curve defining the "hill" corresponding to the peak, the heart beat 
will be properly modelled. The peak searching procedure is illustrated in Fig. 3. 

First the index for the maximum peak is found, Fig. 3 a. Secondly, the area of the 
"hill" is estimated, Fig. 3b, and the complete "hill" is deleted from the spectrum, Fig. 
10 3 c. The procedure is repeated with the modified spectrum, until a maximum number 
of peaks has been found, or the peak level is below a pre-defined threshold. For 
certain heart sounds as many as 50 peaks corresponding to 50 sinusoids may be 
found. 

15 The result from the peak search is a table of peaks for each position of the window 
function. The table corresponding to window position m will be termed spec ra. The 
distance between each window position is T seconds. This is illustrated in Fig. 4. 

The peaks in spec m now have to be matched together with peaks in spec m+1 , in 
20 order to construct a frequency track that the total phase function Q.(t) must follow. 
Many different solutions for this matching problem have been proposed in the 
literature, as no optimal solution exits. The matching strategy used in the present 
invention is based on the following three rules: 

- high energy peaks are matched first 
25 - tracks are not allowed to cross. 

- the frequency differences between two matched peaks must be below certain 
limits. 

First, the procedure finds the maximum peak in spec m. This peak has to be matched 
30 with a peak in spec m+L For this purpose, the procedure searches spec m+1 for the 
maximum peak within lower and upper limits. These limits are controlled by two 
factors. The distance to the nearest peak in spec m and the maximum allowed 
frequency change during T seconds. Before a match is accepted, it must be verified, 
i.e. it must not cross a previously found match. When the final match is verified, the 
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two peaks are no longer visible to the matching procedure. Unmatched peaks in 
spec m are declared dead, and unmatched peaks in spec m+1 are declared bom. The 
procedure is illustrated in Fig. 5, and an example of the matching process is 
illustrated in Fig. 6. 

5 

The synthesis will be designed for T seconds at a time, that is from window position 
m to window position m+1, and this T second reconstructed signal will be termed 
frame m, with a frame length T. The number of born, continued, and dead tracks in 
frame m will be termed Ll(m), L2(m), and L3(m). For each frame ra, the born, 
10 continued and dead tracks are synthesized separately and added to form the final 
signal: 



Eq. (4): 

y (0 = j>i(0+j>2(0 + j> 3 (') 

Ll(m) L2(m) L3(m) 

= £ 4\(0cos(Q™(0)+ X ^S(0oos(n£(0)+ Z < 3 (0cos(n™(0) 

kl = 1 £2 = 1 fc= 1 
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Each track k in frame m has to be synthesised by a sinusoidal function 
,4*(t)cos(Q*<»- Th e track is up to now only specified by the start and end points. This 
is illustrated in Fig. 7. 

20 A m k is the amplitude of the Fourier coefficients, co m * is the frequency of the Fourier 
coefficient and (p m * is the phase offset for the current track. The phase offset is 
increased by an average phase step dependent on the start and stop frequency of the 
track (see Fig. 8). The original model is based on the phase of the Fourier 
coefficients. The amplitude is interpolated by A^t), and the instantaneous phase is 

25 interpolated by Q^O- The interpolation functions are only constrained by the start 
and stop conditions: 



Eq. (5): 





A k 


MT) = 


m + 1 

A k 


n k (0) 


m 

= <Vk 


n k (T) = 




*V(0) 


m 


n k '(T) = 


m+ 1 
= ™k 



The amplitude function is not critical, so a linear interpolation function will be used 
for Ak(t). However, the instantaneous phase function is more critical because 
any discontinuity in the phase progress will deteriorate the perceptual quality of the 
5 final sound. In the following, the synthesis of born, continued, and dead tracks will 
be designed separately. 

In the case of a continued track, a second order polynomial will be used to 
interpolate a smooth phase progress during a frame: 

10 

Eq. (6): 

2 

Cl(t) = a } + a 2 • t+ a 3 • t 



This is illustrated in Fig. 8. 
15 The coefficients are determined by inserting the conditions from Eq. (5) into Eq. (6): 
Eq. (7): 

Q(0) = a x = <p£ 
Q ! (0) = a 2 = <o™ 

m , m + 1 ™ m+1 r* m 

9 m CO/ + 00/ °H — COjl 

C1{T) = a i +a 2 -T+a 3 -T 1 = qff+T- * g * a 3 = 2r 

The linear interpolation of the amplitude function is achieved by: 

20 Eq. (8): 

.m An + 1 

A m k {t) = A m k -t A J^- 



The final equations for the synthesis of a continued track k in frame m is given by: 
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Continued track 
Eq. (9): 

m + 1 m 
* m .m + 1 



5 Two special cases remain to be considered. When a track is declared dead or &ora, 
we do not have both start and end conditions for them. 



In the case of a born track, an instantaneous phase function will be initiated from the 
begining of frame k, with A m k = 0, (o m k = <o m k + \ q> m k = 0, and the second order 
10 instantaneous phase function is reduced to a linear function: 



Born track: 
Eq. (10): 

Qwi/A , m + 1 

k (0 = t ■ <o k 

.m + 1 

When a track is declared dead, the track is terminated at the end of frame m using 
1 5 A m k 1 = 0, ay V 7 = (o m k , and (p m * + 7 = T - <p m k , and the second order instantaneous 
phase function is again reduced to a linear function: 



Dead track: 
Eq. (11): 

A m 

A m k (t) = -f(T-t) 

20 

In order to explain more fully the sinusoidal model according to the invention it is 
prepared for implementation in a mathematical simulation program, such as MatLab 
(TM), and the interpolation functions are hence converted to discrete index. The 
discrete time version is obtained using the following substitutions: 
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Eq. (12): 

t = n • A 
T = N ■ A 

m 

m = Vk # 2TC 

™ k N DFT A 

These substitutions are now inserted into the main equations for the sinusoidal 
model. 



Eq. (13): 

m 4- 1 



m + 1 m 
v _k 271 v k 271 



= 9 - + j£_ . . MA + N DFT' * N DFT* n 2 A 2 

^ ^ Fr A 2(N-A) 

m , /n 271 , 2 , m + l n 

= <p k +n-v k -- + n -(v* - v k)TTJ7 

J\ DFT IV ^DFT 



The time index n is running from 0 to N. The discrete amplitude function is given by: 
10 Eq. (14): 

A k [n] = ^(/i-A) - ^ A- — ^ -/2 

The final equations for a discrete-time continued track is: 



Discrete continued track: 
Eq. (15): 

r i m . m 27Z . 2 , m + \ m N 71 

- q> k +n-v k ■ + n -(v* -^TTU 

A m k \n^Al-n A J^_ 

15 

The discrete function for a bom and dead track can likewise be found following the 
same approach: 
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Discrete born track: 
Eq.(16): 

m + 1 

a k [n] = "A- — •— = n-v k ■ — 

1 DFT a IW DFT 
.m + 1 .m + 1 

<[«] = «A = «-^-- 
* L J N-A N 



Discrete dead track: 
Eq.(17): 

m 

Cl k [n] = n ■ A • — --- + <P k = n ■ v k • r- + cp A 

N DFT A ^ DFT 



The result of coding a heart sound with severe systolic murmurs is illustrated in Figs. 
9, 10, and 1 1 . Fig. 1 1 shows the track pattern that is used to reconstruct the heart 
10 sound and clearly illustrates how tracks are born, are continued, and die, according to 
the spectral content of the original heart sound, and in this respect Fig. 1 1 is a 
practical demonstration of the principle shown in Fig. 6. 

Using the novel sinusoidal model described above, we are now in a position to make 
1 5 modifications to the original signal, which would normally be very difficult to 
obtain. As the model contains a very precise description of how the frequency 
content of the signal is changing over time, it provides a good foundation for making 
signal transformations like time stretching (without changing the frequency 
structure), frequency stretching (without changing the time progress), and joint time 
20 frequency stretching. 

The sinusoidal model may be extended to a time-frequency scaled sinusoidal model 
with small modifications. The joint time frequency scaling is obtained by mapping 
track information at time t and frequency ay in the original model to a new time 
25 position f = p * t, where p is the time scaling factor, and new frequency position 
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co~ = q ■ 1, where q is the frequency scaling factor. A pure frequency scaled 

model, may be obtained by setting p = 1 , and a pure time scaling model may be 
obtain by setting q - 1 . The procedure is illustrated in Fig. 1 1 . 

The time scaling is implemented in the sinusoidal model by changing the frame T to 
5 p * Tin the interpolation functions Q(/) and A(t). Frequency scaling can be 

implemented by shifting the peaks found by the peak searching procedure from co to 
q - co: 



T= pT 
= q • co 

Eq. (18): 

10 When q • co crosses half the sampling frequency, the track is removed. Every time a 
born or continued track has been synthesized, the phase offset <p m * + 1 has to be 
updated to: 



Eq. (19): 



m + 1 _ 0 m /T ^ 



1 5 Other frequency modifications will also be possible. Heterodyne modification is 
obtained by shifting peaks from co to co + r, where r is the amount of the spectrum 
movement. 

The time-scaled sinusoidal model may form the foundation for a method for auto- 
20 scaling a heart sound signal, which scales arbitrary heart sounds to 60 beats per 
minute (bpm). This autoscaling capability is obtained using a scale factor p: 
Eq. (20): 



h 

P= 60 
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where h is the heart rate. Because of the sample rate in the Fourier transforms 
there is no problem in following a changing heart rate dynamically. The usefulness 
lies in the fact that a stressed heart will display a change in the spectral content of the 
heart sound, and focussing on this phenomenon may be simplified if the 
5 simultaneous increase in heart rate may be disregarded. 

In Figs. 13 to 18, the results of applying the above described procedure are shown 
when used on a heart sound. 

10 One particular use of the procedure relates to the transmission of very low frequency 
signals over a channel with a limited bandwidth centred much higher than than the 
low frequency signals, such as heart sound signals to be transmitted over a traditional 
telephone line with the bandwidth 300-3400 Hz. Traditionally this could be obtained 
by heterodyning only, in which harmonic relationships between partials are 

1 5 destroyed. According to the present invention scaling to the frequency range of 
interest will preserve such harmonic relationships. 

The basic parameters that apply to the procedure described above when used with 
arbitrary heart sounds may obviously be chosen by the skilled person according to 
20 specific needs. However, the examples reproduced here are based on the following 
parameters: 



Sampling frequency: 1kHz 

Window length: 64 samples 

25 Window shift: 4 samples 

Window type: Gaussian 

Max. number of sines: 50 

Max. frequency shift of tracks between neighbouring windows: 80 Hz 



30 It will be understood that once the signal has been converted to digital representation 
of data, its manipulation may take place in dedicated processors, RISC processors or 
general purpose computers, the outcome of the manipulation being solely dependent 
on the instructions performed on the data under the control of the program written for 
the processor in order to obtain the function. The physical location of the data at any 
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one instant (i.e. in varying degrees of processing) may or may not be related 

to a particular block in the block diagram, but the representation of the invention in 
the form of interconnected functional blocks provides the skilled person with 
sufficient information to obtain the advantages of the invention. 

5 

The foregoing description of the specific embodiments will so fully reveal the 
general nature of the present invention that others skilled in the art can, by applying 
current knowledge, readily modify or adapt for various applications such specific 
embodiments without undue experimentation and without departing from the generic 

10 concept, and therefore, such adaptations and modifications should and are intended 
to be comprehended within the meaning and range of equivalents of the disclosed 
embodiments. It is to be understood that the phraseology or terminology employed 
herein is for the purpose of description and not of limitation. The means, materials, 
and steps for carrying out various disclosed functions may take a variety of forms 

1 5 without departing from the invention. 

Thus, the expressions "means to ... " and "means for ...", or any method step 
language, as may be found in the specification above and/or in the claims below, 
followed by a functional statement, are intended to define and cover whatever 

20 structural, physical, chemical, or electrical element or structure, or whatever method 
step, which may now or in the future exist which carries out the recited functions, 
whether or not precisely equivalent to the embodiment or embodiments disclosed in 
the specification above, i.e., other means or steps for carrying out the same function 
can be used; and it is intended that such expressions be given their broadest 

25 interpretation. 



